% Design of LP BP y HP shelving filters for offline audio WAV processing.
clear all
close all
clc
fs=44100;
fcl=200;
fcm=100;
fbm=100;
Q=fbm/fcm;
fch=2000;
ts=1/fs;
GL=-15;
GM=0;
GH=0;

[COI,FS,nbits]=wavread('01-Change of ideas.wav');


%Low frecuency Shelving.
KL=tan(pi*fcl/fs);

if GL>0 
    VL=10^(GL/20);
    B0L=(1+KL*(2*VL)^0.5+VL*KL^2)/(1+KL*2^0.5+KL^2);
    B1L=2*(VL*KL^2-1)/(1+KL*2^0.5+KL^2);
    B2L=(1-KL*(2*VL)^0.5+VL*KL^2)/(1+KL*2^0.5+KL^2);
    A1L=2*(KL^2-1)/(1+KL*2^0.5+KL^2);
    A2L=(1-KL*2^0.5+KL^2)/(1+KL*2^0.5+KL^2);
elseif GL<0
    VL=10^(-GL/20);
    B0L=(1+KL*(2)^0.5+KL^2)/(1+KL*(2*VL)^0.5+VL*KL^2);
    B1L=2*(KL^2-1)/(1+KL*(2*VL)^0.5+VL*KL^2);
    B2L=(1-KL*(2)^0.5+KL^2)/(1+KL*(2*VL)^0.5+VL*KL^2);
    A1L=2*(VL*KL^2-1)/(1+KL*(2*VL)^0.5+VL*KL^2);
    A2L=(1-KL*(2*VL)^0.5+VL*KL^2)/(1+KL*(2*VL)^0.5+VL*KL^2);
else
    B0L=1;
    B1L=0;
    B2L=0;
    A1L=0;
    A2L=0;
end

%High frecuency Shelving.
KH=tan(pi*fch/fs);

if GH>0 
    VH=10^(GH/20);
    B0H=(VH+KH*(2*VH)^0.5+KH^2)/(1+KH*2^0.5+KH^2);
    B1H=2*(KH^2-VH)/(1+KH*2^0.5+KH^2);
    B2H=(VH-KH*(2*VH)^0.5+KH^2)/(1+KH*2^0.5+KH^2);
    A1H=2*(KH^2-1)/(1+KH*2^0.5+KH^2);
    A2H=(1-KH*2^0.5+KH^2)/(1+KH*2^0.5+KH^2);
elseif GH<0
    VH=10^(-GH/20);
    B0H=(1+KH*2^0.5+KH^2)/(VH+KH*(2*VH)^0.5+KH^2);
    B1H=2*(KH^2-1)/(VH+KH*(2*VH)^0.5+KH^2);
    B2H=(1-KH*2^0.5+KH^2)/(VH+KH*(2*VH)^0.5+KH^2);
    A1H=2*(KH^2/VH-1)/(1+KH*(2/VH)^0.5+(KH^2)/VH);
    A2H=(1-KH*(2/VH)^0.5+(KH^2)/VH)/(1+KH*(2/VH)^0.5+(KH^2)/VH);
else
    B0H=1;
    B1H=0;
    B2H=0;
    A1H=0;
    A2H=0;
end

%Mid frecuency Shelving.
KM=tan(pi*fcm/fs);

if GM>0 
    VM=10^(GM/20);
    B0M=(1+VM*KM/Q+KM^2)/(1+KM/Q+KM^2);
    B1M=2*(KM^2-1)/(1+KM/Q+KM^2);
    B2M=(1-VM*KM/Q+KM^2)/(1+KM/Q+KM^2);
    A1M=2*(KM^2-1)/(1+KM/Q+KM^2);
    A2M=(1-KM/Q+KM^2)/(1+KM/Q+KM^2);
elseif GM<0
    VM=10^(-GM/20);
    B0M=(1+KM/Q+KM^2)/(1+VM*KM/Q+KM^2);
    B1M=2*(KM^2-1)/(1+VM*KM/Q+KM^2);
    B2M=(1-KM/Q+KM^2)/(1+VM*KM/Q+KM^2);
    A1M=2*(KM^2-1)/(1+VM*KM/Q+KM^2);
    A2M=(1-VM*KM/Q+KM^2)/(1+VM*KM/Q+KM^2);
else
    B0M=1;
    B1M=0;
    B2M=0;
    A1M=0;
    A2M=0;
end
%  z = tf('z',ts);
%  H=((B0L+B1L/z+B2L/z^2)*(B0H+B1H/z+B2H/z^2)*(B0M+B1M/z+B2M/z^2))/((1+A1L/z+A2L/z^2)*(1+A1H/z+A2H/z^2)*(1+A1M/z+A2M/z^2));

%  bode(H)
LPB=[B0L B1L B2L]; LPA=[1 A1L A2L];
BPB=[B0M B1M B2M]; BPA=[1 A1M A2M];
HPB=[B0H B1H B2H]; HPA=[1 A1H A2H];
COIX=filter(LPB,LPA,COI);
COIX=filter(BPB,BPA,COIX);
COIX=filter(HPB,HPA,COIX);
%  figure
%  subplot(2,1,1); plot(COI(1:5*FS));
%  subplot(2,1,2); plot(COIX(1:5*FS));